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We present a method for the calculation of asymptotic operators for non-equilibrium steady-state 
quantum systems. The asymptotic steady-state operator is obtained by averaging the corresponding 
, operator in Heisenberg representation over infinitely long time. Several examples are considered to 

demonstrate the utility of our method. The results obtained within our approach are compared to 
those obtained within the Schwinger-Keldysh non-equilibrium Green's functions. 
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I. INTRODUCTION 



If we place a finite quantum or classical system into the contact with several different macroscopic thermal baths 
or particle reservoirs and wait for some time, which is much longer than typical relaxation time of the system, the 
& ■ system will reach a non-equilibrium time-independent steady state. Like an equilibrium represents stationary state 
^| of a closed system, non-equilibrium steady state is the time-inviriant state of an open system. The steady state is 
established by the delicate balance between irreversible processes and the driving forces produced by the reservoirs. 
Non-equilibrium steady-state systems are ubiquitous and their theoretical description has been a challenging fun- 
damental problem for many years They are also of significant practical interest for various nanotechnological 
and biological applications. Examples of steady-state non-equilibrium systems include quantum contacts Q, molec- 
ular motors [H, molecular junctions and nanowires low dimensional heat conducting quantum and classical 
systems @, H| 0|- Recent experiments on these systems have revealed a wealth of interesting new non-equilibrium 
^ | phenomena such as non-diffusive heat transfer, negative differential resistance, stochastic switching and hysteresis of 
electric current, and direct harnessing of thermal fluctuations [ToL ITU H3|. 

The understanding of the fundamental mechanisms as well as the interpretation of these experimental observations 
require the development of new theoretical methods for the description of non-equilibrium quantum systems. Many 
' theoretical approaches have been developed to deal with non-equilibrium steady-state systems, such as Keldysh- 
Schwinger Green's functions [Tj3, 0, H, EH, E3, 0, El , scattering theory based approaches [20L l2lj | , different variations 
\SD • of Zubarev method of non-equilibrium statistical operator [U, IH, HH, Hf| |2q |. and constrained current methods 
S3, EE H, US HI, HI H H|. Most of the methods are based on calculations of non-equilibrium averages of the 
\Q • operators of interest. Although average values, such as correlation functions, densities, currents, are very important 
and directly related to the measurable quantities, it is very useful to have explicit expressions for asymptotic operators. 
ON . Here we propose a general theoretical approach which directly defines and computes asymptotic non-equilibrium 
steady-state operators. 



II. ASYMPTOTIC OPERATORS 

Let us consider a quantum mechanical observable which corresponds to a Hermitian operator A. We assume the 
existence of a unique steady state for our system. It implies the existence of the limit t — > Too for observable 
quantities. It means that the asymptotic value of the operator A in Heisenberg representation is uniquely defined and 
can be computed as 

A= lim e mt Ae- lHt , (1) 

t — >+oo 

where the Hamiltonian H includes the system Hamiltonian Hs, the bath Hamiltonian Hb, and the system-bath 
interaction Hsb- 



H = Hs + Hb + Hsb ■ 



(2) 
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We denote all thermal baths and particle reservoirs collectively by B. The Hamiltonian H is assumed to be time- 
independent. We begin calculation of the limit |T]) by the use of the following self-evident identity, which is valid for 
any time-dependent operator 0(t): 



<>in = ow fdt'^p-. (3) 



Applying this identity to eq.([T]), we get 



— f + °° d 

A = A+ dt—e lHt Ae~ lHt . (4) 
jo dt 

The existence of the unique steady state implies the convergence of the integral in eq. (U]) , therefore we can introduce 
the infinitesimal factor ?y: 

_ r + °° d 

A = A+ lim / die-"'— e lHt Ae~ lHt . (5) 
v^+QJo dt 

This regularization of the integral is equivalent to the assumption that there exist some dissipation mechanisms in 
our system which we do not specify explicitly but which are efficient enough to lead to the unique steady state. The 
integration by parts of eq. (O yields 

r+oo 

A = lim n / dte-^e lHt Ae- lHt . (6) 
n^+o 'Jo 

This expression for the non-equilibrium asymptotic operator is interesting physically since the integral is multiplied 
by the infinitesimal number rj. Therefore, only singular terms proportional to l/rj give non- vanishing contribution 
from the integral to A. Using Abel's theorem [22] 

lim- / dtf(t)=]imri dte~*f(t), (7) 

we see that eq.© yields the part of the operator which is conserved after averaging over the infinitely long time. 
Eq.® is the definition of the asymptotic non-equilibrium steady-state operator. It will be the starting point for our 
calculations. The definition ([6]) has been initially introduced by Kubo in his work on irreversible thermodynamics 
(so-called time invariant part of the operator) [35| and Zubarev during his development of non-equilibrium statistical 
operator method [13] ■ It was used Grandy (36| and Hershfield [23{ in their formulations of steady state non-equlibrium 
statistical mechanics. Recently, Tasaki and Takahashi also employed asymptotic operators ^ to continue the advances 
of Zubarev method in application to transport in quantum junctions [241 ] . The time invariant current defined via ([7]) 
was used by Bokes, Mera, and Godby as a term constrained by the Lagrange multiplier in their development of 
variational transport theory[31]. In this paper, we develop a general practical method for computing asymptotic 
non-equilibrium operators. 

To carry out our program it is convenient to represent the Hamiltonian and operator A in the second quantization 
form. We assume that the Hamiltonian and the operator A are quadratic in creation and annihilation operators: 

H = Y J H l] a\a ]1 (8) 



A = J2 A iM a r ( 9 ) 

ij 

The creation at and annihilation a operators obey the standard commutation (for bosons) and anticommutation (for 

fermions) relations. The methods developed in this paper are applicable to both, Fermi and Bose, systems. With the 
use of the following operator identity 

e mt Ae -mt = J2 ( e ™Ae- lHf ) y . a\a h (10) 

ij 
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where H and A are matrices with matrix elements iJy J8]) and Ay ([9]), the expression for the asymptotic operator 
© becomes 



r+oo 

A= lim^y rft^] | e -<(«-H-ii)t Ae «( u -H+i,) 



(11) 



This integral can further be transformed by inserting the delta function 6(t — t r ) = /_ °° duie^^ * ) and integrating 
over t and t': 



A= lim -L 

7T 



+ oo 



(12) 



where G(o>) = [cj — H + i?/] -1 is the total Green's function. This is the general formula for the asymptotic steady-state 
non-equilibrium operator, which is applicable to any systems with quadratic Hamiltonians. 

The indexes i and j run over the system s and bath b single particle states in cqs.® and ([9]). Therefore the total 
Hamiltonian matrix can be written in the following form 



H 



(13) 



where H5 is the matrix of the system Hamiltonian, H# is the matrix of the bath Hamiltonians, and Hsb is the 
matrix of the system-bath interaction. Likewise, we can partition the matrix A 



A B A SB 
A B s A s 



and the Green's function 



G(w) = [u - H + w?]" 1 



Gb(u) Gsb(w) 
Gbs(w) Gs(w) 



(14) 



(15) 



The block matrix elements of the Green's function, which can be computed by the Frobenius formula. [37j| have the 
following form 



GsM = [(w + irj)I s - H s - E B (c< 



(16) 



GsbM = -G s (w)H SB g B (w), 



(17) 



Gb(w) = gsM + g B (w)H BS Gs(w)H lSB g B ((x)). 
Here we introduced the Green's functions for the baths 



and the corresponding self-energies 



gs(w) = [(cj + w?)I B - H B ] 



£b(w) = H S BgB(w)H BS . 



(18) 



(19) 



(20) 



Note that the baths are assumed to be macroscopically large, so they have continuum energy levels. Therefore, when 
T) tends to zero, the Green's function of the baths gs(^) has a pole on the real energy axis. The Green's function 
Gs(w) is not singular in general, although special care should be taken to deal with the bound states [Hj]. 

We consider two practically important choices of the operator A at time t = 0. First, A is defined in the system 
space only 



A s = lim i / GL,V[G>)A s G(u;)]..a a, = lim i / dw V [G^(w)A s G SB (. 



w)] 66 ,oJat/. (21) 
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Second, A mixes the Fock spaces of the system and baths: 

p+oo 



/-t-oc /'-hoc 
du V [G*(w)A SB G(w)] y ajoj = lim J / du; V [G£ s MA SB G B sML, a£<v 

/ + OG Z' + OO 

dw5^[GSHAsflG B (w)] H) ,atay+ hm/ / d w ^ [G BS HA SB G B H] bb , 4^22) 



66' 66' 



Here the indexes b and s refer to the single-particle states of the baths and the system, respectively. We have retained 
only those terms under the integrals which are ~ l/((w — Eb) 2 + rf). Less singular terms give zero contributions 
upon the integration and multiplication by rj. With an eye on the averaging of asymptotic operators over the initial 
equilibrium density matrix, which is the product of equilibrium system and bath density matrices, the expression for 
Asb can further be simplified: 

A ss «lim// d w ^[G BS HA SB G B H] bb ,a^. (23) 



III. COMPARISONS WITH NON-EQUILIBRIUM GREEN'S FUNCTION METHOD 

In order to compare our approach with existing theoretical methods, we consider a finite quantum system with 
discrete energy levels and non- interacting fermions (quantum dot, atom, atomic cluster, molecule, or atomic Fermi-gas 
in a trap) which is attached to two macroscopic reservoirs. Such kind of systems can be described by the tunneling 
Hamiltonian [391 ]. which has the form 

H = H Q + T. (24) 
Here Hq is the Hamiltonian of the isolated reservoirs and the system, 

H Q = £ 6 a fe a 6 +y^g a <4 a s- ( 25 ) 

b=l,r s 

The indexes I and r refer to the continuum single-particle states in the left and right reservoirs, respectively, and 
s refers to the discrete single-particle states of the system. The tunneling interaction couples the reservoirs to the 
system, 

T = y^t6 s (a|q s + a\ab)- (26) 

6=Z,r s 

We deal with fermions in this example, therefore the creation a* and annihilation a operators obey the standard 
anticommutation relations. We assume that the system is isolated from the reservoirs at t < . Therefore the system 
is in equilibrium at t < and it is described by the density matrix 

Here Nl (Nr) are the number of particles operators for the left(right) reservoirs and \il (vr) are the corresponding- 
chemical potentials. Then at t = we turn on the reservoir-system coupling T. 
First, we compute the asymptotic single-particle density matrix for the system 



The use of formula (|21|) gives 



lim e lHt a\,a s e- lHt . (28) 



/-f oo 
dw£ [G BS HL, [GsbML, "I'v, (29) 



where b and b' run over all single-particle states of the left and right reservoirs. Substituting the explicit expressions 
for the Green's functions (fT7|) . we get 



r/—>+Q 7T 



lim// ^^[g^MH ss GSM] bs ,[G s HH SB g B H] sb ,a \a v . (30) 
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Averaging the asymptotic operator (|30[) over the equilibrium initial density matrix, we obtain the following expression 
for the steady-state single-particle density matrix: 



1 r+°° 

Tr{p n ss ,}=— dco J2 /^-Mb)[G 5 M(S| j H-S b H)G^M] ss ,, 

•J —OO d r d 



(31) 



B=L,R 



f(u>) being the Fermi-Dirac occupation numbers. The same equation is obtained by the method of non-equilibrium 
Green's functions through the lesser Green's function G < [io| . 

Now we calculate the asymptotic current operator. The current operator is defined by the continuity equation: 



J = j t N L = i[H,N L ] =iY,Ti s {a\a s - a\ ai ), 



(32) 



Nl being the particle number operator for the left reservoir [4l| . 
We define the asymptotic value of the current operator as 

7= lim e lHt Je 

t — >+oc 



iHt 



The use of eq. ([23|) leads to 



J = lim ■ 

>)^+0 7T /_. 



+oo 



b av, 



(33) 



(34) 



where the matrix, which represents the current, is 



1,1/ 



T LN 
-T NL 




(35) 



Substituting the expressions for Gb(oj) and Gss(w) into eq. (|34|) . and averaging over the equilibrium initial density 
matrix (|2T[) results into the following expression for the steady-state current: 



2tt 



duj{n(uj - fx L ) - n(u - ^fl)) Tr [( s l( w ) ~ s z,( w )) G iVAr( w )( s fl( w ) ~ s /j(^))Gatjv(w)] 



(36) 



This is the so-called Landauer formula, which is typically obtained via the formalism of non-equilibrium Green's 
functions HO. 



IV. CONCLUSIONS 



We define a non-equilibrium steady state as an asymptotic state of a finite quantum system which is connected 
to several different macroscopic thermal baths or particle reservoirs. We present a general method to compute 
asymptotic steady-state operators for such systems. The asymptotic operators can be computed exactly for the 
systems with quadratic Hamiltonians. For the tunneling Hamiltonian we recover the standard results obtained by the 
non-equilibrium Green's function method, thereby demonstrating equivalence of the two approaches. For example, 
we derive Landauer formula by a direct calculation of the asymptotic steady-state non-equilibrium current operator. 
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